library(dpylr)
library(monocle)
library(Seurat)

Oligos <- readRDS("/media/annawilliams/Data/SarahDatastore/Nina/HumanOL_Oligos")

#Extract data, phenotype data, and feature data from the SeuratObject
oligo.data <- as(as.matrix(cultured.oligos@assays$RNA@data), 'sparseMatrix')

pd <- new('AnnotatedDataFrame', data = cultured.oligos@meta.data)

fData <- data.frame(gene.short.name = row.names(oligo.data), row.names = row.names((oligo.data)))
fd <- new ('AnnotatedDataFrame', data = fData)

#Construct monocle cds
cultured.oligos.cds<- newCellDataSet(oligo.data,
                                     phenoData = pd,
                                     featureData = fd,
                                     lowerDetectionLimit = 0.5,
                                     expressionFamily = negbinomial.size())
# this opens a new wqindow with a table with all information
View(pData(cultured.oligos.cds)) 



cultured.oligos.cds <- estimateSizeFactors(cultured.oligos.cds)
cultured.oligos.cds <- estimateDispersions(cultured.oligos.cds)

#we can do the pseudotime based on differentially expressed genes between the clusters
#I have used once 1000 genes and then 1500 but it looks very similar, so no real difference, also gives you 5 states

clustering_Oligos <- differentialGeneTest(cultured.oligos.cds,
                                          fullModelFormulaStr = '~seurat_clusters',
                                          cores = 15, verbose = T)

# select genes that are significant at an FDR >10%
genes_cultured.oligos.cds <- subset(clustering_Oligos, qval < 0.1) 


# my_ordering_genes <- row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)][1:1000]
cultured.oligos.cds <- setOrderingFilter(cultured.oligos.cds, ordering_genes = genes_cultured.oligos.cds)
cultured.oligos.cds <- reduceDimension(cultured.oligos.cds, method = 'DDRTree', cores=6, verbose = TRUE)
cultured.oligos.cds <- orderCells(cultured.oligos.cds)

plot_cell_trajectory(cultured.oligos.cds, color_by = "RNA_snn_res.0.8")
plot_cell_trajectory(cultured.oligos.cds, color_by = "RNA_snn_res.0.8") + facet_wrap(~RNA_snn_res.0.8)

plot_cell_trajectory(cultured.oligos.cds, color_by = "State") + facet_wrap(~State)
plot_cell_trajectory(cultured.oligos.cds, color_by = "Pseudotime")